%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This Matlab program is for Figure 1 from the Supplemental Appendix of 'The explicit formula for the Hodrick-Prescott filter in finite sample' by Adriana Cornea-Madeira
% HP filter weights computed using our formulae from Corollaries 6 and 7, 
% for n even
clear all; clc;
alpha = 1600; l=2;
nnn = (100:100:10000)';
Kk = 50; % number of samples
t = zeros(length(nnn),1);
for in=1:length(nnn)
tic
n = nnn(in,1); % sample size
m = n/2;
fun = @(x,z) x.*z;
i = 1:m;
x = (-1).^(i-1);
z = ones(m,1);
J = bsxfun(fun,x,z);

fun = @(i,j) sqrt(2/(n+1))*sin(i.*j);
j = i'*pi/(n+1);
T1 = bsxfun(fun,i,j);

fun = @(i,j) sqrt(2/(n+1))*sin(i.*j*2);
M2 = bsxfun(fun,i,j);

i = i';
j = i;

a1 = sin((2*j)*pi/(n+1));
b1 = sin(2*(2*j)*pi/(n+1));
c1 = 1+alpha*(2-2*cos((2*j)*pi/(n+1))).^2;

fun = @(a,b) (2/(n+1))*a.*b;
a = (2*a1-b1)./c1;
b = a';
G2 = bsxfun(fun,a,b);

k2 = 2*alpha/(1-2*alpha*sum((2/(n+1))*(2*a1-b1).^2./c1));
E = k2*M2*G2*M2';

% eigenvalues
b = cos(j*pi/(n+1));
gam1 = 1./(1+4*alpha*(1-b).^2);
gam2 = 1./(1+4*alpha*(1+flipud(b)).^2);

X = T1.*J; 

V1 = T1*bsxfun(@times,gam1,T1) + X'*bsxfun(@times,flipud(gam2),X);
V2 = T1*bsxfun(@times,gam1,fliplr(X'))+(-1)^l*X'*bsxfun(@times,flipud(gam2),fliplr(X')).*J;
j = i-1;
a2 = sin((2*j+1)*pi/(n+1));
b2 = sin(2*(2*j+1)*pi/(n+1));
c2 = 1+alpha*(2-2*cos((2*j+1)*pi/(n+1))).^2;
k1 = 2*alpha/(1-2*alpha*sum((2/(n+1))*(2*a2-b2).^2./c2));
    
fun = @(i,jj) sqrt(2/(n+1))*sin(i.*jj);
jj = (2*j+1)'*pi/(n+1);
M1 = bsxfun(fun,i,jj);
fun = @(a,b) (2/(n+1))*a.*b;
a = (2*a2-b2)./c2;
b = a';
G1 = bsxfun(fun,a,b);
  
D = k1*M1*G1*M1';
     
H = V1+D+E;
DE = D-E;
    
% the inverse
Finv_new = [H V2+fliplr(DE);rot90(V2,2)+flipud(DE') rot90(H,2)]; 
for k=1:Kk
  % local linear trend model (LLT)
  beta = zeros(n,1);
  tau = zeros(n,1);
  y = zeros(n,1);
  sigma_c = 2;
  sigma_eta = 0; %smooth trend model
  sigma_zeta = 0.025*sigma_c; 
  beta(1) = randn(1,1); tau(1) = randn(1,1); y(1) = randn(1,1);
  epsi_1 = randn(n,1); epsi_2 = randn(n,1); epsi_3 = randn(n,1);
  for i = 2:n
      beta(i) = beta(i-1) + sigma_zeta*epsi_1(i,1);
      tau(i) = tau(i-1) + beta(i-1) + sigma_eta*epsi_2(i,1);
      y(i) = tau(i) + sigma_c*epsi_3(i,1);
  end
tau_est = y - Finv_new*y;  
end
t(in)=toc;
end
figure(1)
plot(nnn,t,'m'); hold on;